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Some efficient Nordsieck integration 
methods for IVPs 


N. Barghi Oskouie, A. Abdi* and G. Hojjati 


Abstract 


In this paper, in continuation of the construction of efficient numerical 
methods for stiff [VPs, we construct type two Nordsieck second derivative 
general linear methods with order p = s, where s is the number of internal 
stages, and stage order gq = p. Implementation of the constructed methods 
with fixed and variable stepsize is discussed which verifies their efficiency. 


Keywords: Stiff differential equations; Nordsieck second derivative general 
linear methods; A- and L-stability; Variable stepsize implementation. 


1 Introduction 


Second derivative general linear methods (SGLMs) for the numerical solu- 
tion of autonomous ordinary differential equations (ODEs) with initial value 
problem 


(1) 
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where f : R™ — R”™ and m is the dimensionality of the system, have 
been studied in recent years. These methods were introduced by Butcher 
and Hojjati in [1] and investigated more in [246] by Abdi and Hojjati. In 
the construction of SGLMs which are extension of general linear methods 
(GLMs) [10,((3H15,((7]), there are a lot of free parameters which allow us 
to construct high order methods with a small number of internal stages to 
reduce computational cost [ll]. 


We recall that SGLMs are characterized by four integers, (p, q,7r,s) where 
p and q are, respectively, order and stage order of the method, r is the number 
of input and output approximations, and s is the number of internal stages. 
Let Yl) = [y/"]s_, be an approximation of stage order q to the vector 
y(@n—1 + ch) = [y(@n—1 + cih)]$_1, and let the vectors FV") = [FY Ve, 
and g(Yl"l) = (o(v")] ’_, denote the stage first and second derivative values, 
where g(-) = f’(-)f(-), respectively. If r = p+ 1, we can assume the input 
and output vectors at the step number n, y!"~!) and yl"), are approximations 
of order p to the Nordsieck vectors 


y(an—1) y(n) 

hy'(#n—1) hy' (xn) 
and 

hPy®) (tn-1) hPy) (xn) 


respectively. In an SGLM used for the numerical solution of (I), these values 
are related by 


Yl = h(A @ Im) f(¥"™) + h?(A @ Im)g(¥™) + (U ® Im)y-4, 
(2) 
yl = h(B@ Im) f(VM) + h?(B® In)g(¥™) + (V @ Im)y4, 


where n = 1,2,...,N, Nh = %—2o, h is the stepsize, and ® is the Kronecker 
product of two matrices. Here A, A € R***, U € R**", B,B € R™®, and 
V € R'™". The coefficients matrix V in the Nordsieck SGLM (@) has the 
form 


v= [vp vo +++ ea 5 V € RO-YDx(r-1)_ In this paper, the methods will 
be restricted to the case where p = q = r — 1 = s and the eigenvalues of V 
are zeros. The latter condition ensures zero-stability of the method. 


An SGLM in Nordsieck form has order p and stage order q = p if and 
only if [Xj 
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U =C— ACK — ACK?, 
V = E-— BCK — BCK?, 


where the matrices C € R°X®@+)), k € R@+Yx@+) and E € R@t)x(e+1) 
are defined by 
2 cP 


Cc Cc 
C:=)}1 Tl a a8 pl 5 K := (0 ey €2 eee €p |, 


and 


E := exp(K) = 001 ae : ; 


i 
1! 
1 


with e; as the jth vector of canonical basis in R?*', respectively. 


The stability behavior of SGLMs is defined using the standard test prob- 
lem of Dahlquist y’ = €y, where € is a complex number. If method () is 
applied to this problem, then the stability matrix is 


M(z)=V+(2B+2°B)(I-zA- 2A) 'U, 


where z = h€ and the stability function of the method is defined as the 
characteristic polynomial of M(z); that is, 


p(w, z) = det (wI — M(z)). 


If p(w, z) = w’!(w — R(z)), the method is said to possess “Runge-Kutta 
stability (RKS)”. GLMs and SGLMs with RKS property were studied by 
Butcher and Jackiewicz in [9, T2414] and by Abdi and Hojjati in BH, re- 
spectively. 


For economical implementation, it is assumed that the matrices A and 
A have a lower triangular form with the same diagonal entries \ and p, 
respectively. SGLMs are divided into four types, depending on the nature 
of the differential system to be solved and the computer architecture that is 
used to implement these methods. Types 1 and 2 are those with arbitrary 
aij and a;;, where \ = up = 0 and A > 0, ps < 0, respectively. Such methods 
are appropriate, respectively, for nonstiff and stiff differential systems in a 
sequential computing environment. Requiring a;; = aj; = 0, cases AX = uw = 0 
and A > 0, uw < 0 lead, respectively, to types 3 and 4 methods which can 
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be useful, respectively, for nonstiff and stiff systems in a parallel computing 
environment. 

The construction and implementation of type 2 SGLMs with p = s+1 
were discussed in [A]. Also, the construction of parallel Nordsieck SGLMs 
with p = s and their order barriers were studied in [5]. These order barriers 
were obtained under the assumption of RKS property. In this paper, we are 
going to construct type 2 Nordsieck SGLMs with p = s, which are efficient 
methods for stiff systems. Also, efficiency of the constructed methods are 
shown by their implementation in a variable stepsize environment. 

Next sections of this paper are organized as follows: In section J, we con- 
struct A- and L[-stable SGLMs in the Nordsieck form with RKS property of 
orders 2, 3, and 4. Considering Nordsieck SGLMs in the variable stepsize 
mode, implementation issues including local error estimation, and stepsize 
control are discussed in section B. Finally, in section M, some results of nu- 
merical experiments on some stiff test systems are presented and compared 
with those obtained by Nordsieck SGLMs of the same order. 


2 Construction of type 2 methods 


In this section, the construction of type 2 Nordsieck SGLMs of order p and 
stage order g = p with some desired stability properties is explained. 


2.1 Methods with p=q=s=r—-1=2 


We construct methods with p = q = s = r—1= 2 and RKS property. We 
look for methods which their stability function has the form [5] 


_ L4+nyztnez? + 132° 


#2) (1—Az— p22)?” 


where 2 
1 + So ngat = exp(z)(1 — Az — 2”)? — C2? + O(2*) 
k=1 
with C as the error constant of the method. For this method to be A-stable, 


it is necessary and sufficient that \ > 0,  < 0, and so that the E(y) is 
non-negative for all real y, where the E-polynomial is defined by 


E(y) = |1— diy + py? |* — | + rity — nay? — gig”. 


By choosing C = 107‘, a detailed calculation shows that 
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E(y) = y*(Eo + Eiy? + Eoy*), 


where 
1247 4997 
Eo = Poy +? = 2 N+ 4p, 
cou Y~ 7500 e 
24970009 19997 4997 4997 
E, = —-——_—_ — 4 + 23 d? — Ay? 
1 ~ 900000000 - 15000 H+ Ze00"T T5000 
34997 
tA? — 2)? 1? — AP u 4 24 8\2 4 — 
Ll LL Ut BAY + BA" pL 7500 Me 
Ey = 7". 


Pairs of (A, #4) with values in domain [0,2] x [—2,0] giving L-stability are 
shown in Figure I. 


0 
0.25 
[Lb 
0.5 
—0.75 
0 0.5 1 1.5 2 
my 


Figure 1: L-stable choices of (A, 4) for p = s = 2 corresponding to C = 10~¢. 


We select a single example, characterized by \ = - w= —2 ande = 


5 
[5 1]. The coefficients of the method are 


r 4 i —~38 _ 3 

5 0 5 0/1 10 40 

967. 4| 494 1), AT17 253 

18750 5| 3375 5 18750 12500 

A| A| U 967 494 1), Az17 253 

= 18750 3375 5 18750 12500 
B\B\|V 


4 
5 

0 1) 0 Oj|O 0O 0 
0} 0 1/40 O 0 
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2.2 Methods with p=q=s=r—-1=3 


We construct methods with p = q = s = r—1=3 and RKS property. We 
look for methods which their stability function has the form [5] 


( ) Lltniztnez? +n323 +naz*+n5z° 
R A— 
(1 — Az — puz?)3 


where 


1+ S- n,z” = exp(z)(1 — Az — uz”)® — C,2* — Coz + O(2°). 
k=1 


Here C, is the error constant of the method and C2 is an arbitrary number. 
The E-polynomial has the form 


E(y) = y®° (Eo + Eiy? + Eoy* + Esy’). 


For these methods to be A-stable, it is necessary and sufficient that ’ > 0, 
pe < 0, and so that Ey + Eyx + Eoa? + E32? + Eya* is non-negative for 
all positive real numbers x, where Eo, £,, E2, E3, and E, are complicated 
expressions in terms of \ and pz. By choosing C; = C2 = 10~+, pairs of (A, 1) 
with values in domain [0,2] x [—2,0] giving L-stability are shown in Figure 
2. 


—0.1 


—0.2 


—0.3 


0 0.5 X 1 1.5 


Figure 2: L-stable choices of (A, j1) for p = s = 3 corresponding to C1 = C2 = 10-4. 


Here, we represent an example, characterized by 


1 1 1. 2 


c= F 3 14 fae 
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The coefficients of the method are 


T 0.5000000000000000 0 0 
A = | 1.4279081052775164 — 0.5000000000000000 0 
1.0000000000000000 —0.3168631901664915 0.5000000000000000 
| —0.0666666666666667 0 0 
Z = | —0.3067166674763493 —0.0666666666666667 0 
—0.0602082721233515  0.0288951398441268 00sec 
1.0000000000000000 +—0.3168631901664915 — 0.5000000000000000 
0 0 1 
ae: 0 0 0 } 
84.1340111524194390 —15.9895442199910120 -37.9511398307057703 | 
[ —0.0602082721233515  0.0288951398441268 —0.0666666666666667 
_ 0 0 0 | 
a= 0 0 1 ? 
0 —1.7866934603873189 20.0458159571414841 


[ 1.0000000000000000 —0.1666666666666667 —0.0444444444444444 0.0006172839506173 
U = | 1-0000000000000000 —1.2612414386108497 —0.2136971453939340 0.0056267104705261 
1.0000000000000000 —0.1831368098335086 —0.0241114076097811—0.0010021824846360 


[ 1.0000000000000000 —0.1831368098335086 —0.0241114076097811—0.0010021824846360 ] 


0 0 0 0 
ae 0 0 0 0 
0 —30.1933336017226565 2.3070365964725901 0 


2.3 Methods with p=q=s=r—-1=4 


In this part, we construct methods of the Nordsieck SGLMs of type 2 with 
p=q=s=r-—1=4andc= (0 0 0 1)’. Setting some free parameters 
in order to make calculation easier, RKS conditions make the coefficients 
matrices of the method to take the following forms 
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+00 0|-4 0 0 j1- 400 
330 0/-4-4 0 Of1-1 400 
5 15 0|-% 1 -§ O [1-2 -2 00 
s11$/-} 1 -1 -#$/1 0 300 
s-115/-; 1 -1 -gll 0 § 00 
0000} 0 0 0 1/0 0 000 
600-6} 2 0 0 4/0 0 000 
1200-12} 7-1 0 6/0 0 000 


_ a 


R(z) = _ 


DLW Jdle 
me 
wo 


consequently, by the Ehle conjecture [19], the method is A-stable. The error 
constant of the method is C = =. 


3 Implementation of the methods in a variable stepsize 
mode 


In this section, we recall some implementation strategies given in [A], to apply 
the constructed methods in variable stepsize environment. 
A Nordsieck SGLM in the variable stepsize mode takes the form 


YM =h,(A@ Im) f(V™!) + h2(A@ In)g(¥!) + (UD(Sn) ® Im)y”4, 


yl] =h,(B® Lape) Tr ho (B ® Lg) + (VD(bn) ® ing 4 
(3) 
where 
D(6n) := diag(1,6n,62,...,62), 6, = hn/hn-1. 


Due to the structure of the matrices V and D(d,,), the matrix VD(6,) has 
one eigenvalue with magnitude one and a zero eigenvalue with multiplicity p, 
for any value of dy. 

To achieve a suitable choice of stepsize for the next step, we first need 
to estimate the leading term in the local truncation error. To do this, we 
approximate h?+1y(?+1)(z,,); so that LTE(an) = Cph?tty®t)(2,) can be 
calculated as an approximation to the local truncation error. 

For the methods with p = s = 2 and abscissas c = [5 1]", we use the linear 
combination of the form 
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est(gn) = Cy (ahf(%) + aghf(Y2) + Bh2g(%1)), 


with 
ay = —8, ag= 8, B = —4, 


For the methods with p = s = 3 and abscissas c = [} 3 1)", we use the 


linear combination of the form 
est(tn) = Cp (ahs) + aghf (Ys) + aghf(Y3) + Bh?g(%1)), 


mn 243 81 
ON) = = ag = —162, a= G9 B = 27. 


For the methods with p = s = 4 and abscissas c = [0 0 0 1]", we use the 
linear combination of the form 


est(tn) = Cp (aihs(¥%) + aghf(Y4) + Bih*g(¥3) + Boh? g(Ya) + wt), 


with 
a, = 72, ag=—72, 8, =48, Bo =24, y=12. 


To control the stepsize, we use the following strategy 
est (tn) < Rtol - max{|lyn|, ||Yn+il|} + Atol, (4) 


where Atol and Rtol are given absolute and relative tolerances, respectively. 
If the control (@) is not satisfied, the current step is repeated with a halved 
stepsize. Otherwise, the current step is accepted and we carry out the next 
step with the new stepsize defined as 


An41 = On+ihn, 


where 


On4t1 = min{ facmac, (, gp OE: ) a : 


lest (n.) |loo 


Here, facmaz and fac are safety factors built into a code to prevent the 
step from increasing too rapidly and to avoid an excessive number of rejected 
steps. We have chosen Atol = Rtol = tol, facmax = 2, and fac = 0.9. 


4 Numerical result 


In this section we present the results of numerical experiments to show effi- 
ciency of the constructed methods in section for fixed and variable stepsize 
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mode using the provided techniques in section B. We consider the following 
test problems: 


e Problem 1. The nonlinear stiff system of ODEs 
yi (x) = —1002y: (x) + 100043 (x), 
Yo(©) = yr (x) — y2(x)(1 + yo(x)), 


whose exact solution is [y:(x), yo(a)|7 = [exp(—2z), exp(—x)]7 and 
x € [0,2]. This problem is stiff with an approximate stiffness ratio of 
10° near to = 0. 


Problem 2. A system of differential equation, is called CUSP, resulting 
from discretization of the diffusion terms by the method of line the 
periodic boundary-value problem [6], is as below 


2 
Sy — 13 + ay+b)+o54, 
G2 — $40.07 +0£3, 
9a — (1 — a?)b—a — 0.4y + 0.035v +0 25, 
where i: 
= = (y —0.7)(y — 1.3). 
vege UH WH Ong—13) 


This problem takes the form 
Yi = —E 1 (y? + aays + b:) + D(yi-1 — 2yi + Yi), 
a; = b; + 0.071; + D(aj_1 — 2a; + ai41), 7=1,2,...,N, 
b; = (1 — a?)b; — a; — 0.4y; + 0.03514 + D(bi_1 — 2b; + bi41), 


us = (yi — 0.7)(y: — 1.8), D=oN’, 
with periodic boundary condition 


Yo (= YN; do = an, bo := by, 


YN41°= 41, ON41:= 41, bn 41 := by. 


We take o = Tre € = 107-4, N = 32, and the initial values as 
(0) =0, a,(0)=—2c0s(2), a,(0) = 2sin(™), 4 =1,2,...,.N 
Yi =U, Ay => COs N 5 ay = 2458S1n WN 5 t= 1,4,..., 5 


with tout = 1.1. 
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e Problem 3. For the ODE case, the Ring modulator problem originates 
from electrical circuit analysis is of the form [13] 


y, = C*(ys — 0.5y10 + 0-5y11 + yra — Ro*y1), 
yg = C7! (yo — 0.5y12 + 0.5413 + yis — R7*y2), 
y3 = Cr ‘(no — q(Up1) ) +q(Ups)), 

Ya = cr1(- ya + q(Up2) — q(Ups)); 

ys = Co *(u2 OU pi) = q(Ups)), 

Yo = cs1(- yis — q(Up2) + q(Ups)), 

yh = Cy "(Ro yr + q(Up1) + (U2) ~ q(Ups) — @Us)). 
¥g = =1y, ti; 

yg = —Li ye, 

Yio = Lye (0.541 — ys — Rg2yio), 

Yi1 = Ly (—0.5y1 + ys — Ryayit), 

Yio = Leg (0.5y2 — ys — Rgzyiz), 

¥i3 = Leg (—0.5y2 + ye — Rgsyis), 

Via = La (- y1 + Uini(t) — (Ri + Ry)y); 


Vis = La (- yo — (Re + Rei)is): 


The auxiliary functions Up,, Up2, Ups, Upa, gq, Vini, and Uing are given 


by 


Up1 = ys — Ys — ¥7 — Vina(t), 
Up2 = —y4 + y6 — y7 — Vina (t), 
Ups = ya t+ ys + y7 + Vina(t), 
Upa = —y3 — Yo + ¥7 + Uinalt), 


q(U) > (eu > 1), 
Uini(t) = 0.5 sin(2000zt), 
Uing(t) = 2sin(20000rt). 


The values of the parameters are 
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Calpxcit®, C,=2x1 €,=10-", ta = 445, 


Le1 = 0.002, Leg =5x10-*, Lo3 =5 x 107+, 
R = 25000, Ry = 50, RaS 363, Rye = 173, 
Rg3 = 50, R, = 50, R. = 600, 


5 = 17.74933332, y = 40.67286402 x 107°, 


with the initial vector yo = (0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)" and t € 
(0, 10-3]. 


4.1 Fixed stepsize experiments 


We first present fixed stepsize numerical results in order to show accuracy 
of constructed Nordsieck SGLMs and validate the order of these methods in 
the integration of stiff differential systems. To do this, we have applied the 
methods of order 2, 3, and 4 to the Problem 1. We have implemented the 
methods with a fixed stepsize h = 1/2", with several integer values of k. The 
results of numerical experiments for type 2 Nordsieck SGLMs are shown in 
the Tables 1], 2, and B. In these tables, we have listed the norm of error 
|e, (Z)|| at the endpoint of integration ¥ = 2 and numerical estimate to the 
order of convergence, p, computed by the formula 


_ log(llen@)II/Mlens2@)I1) 
log(2) ; 


where e;,(%) and ep; 2(%) are errors corresponding to stepsizes h and h/2 for 
Nordsieck SGLMs. Also, to show that the constructed methods are com- 
petitive with the efficient existing methods, we have reported the numerical 
results of type 2 methods which have been constructed in [2]. 


Table 1: Numerical results for Nordsieck SGLMs of order p = q = 2. 


Type 2 method of order 2 [2] Nordsieck SGLM of order 2 


" en(Z) Pp en(Z) Pp 
10. «1.78 x 1071! 1.35 x 10718 

11 5.44 107!2 1.71 3.31 x 1071! 1.88 
12 1.43 x 107-1? 1.93 8.18 x 10712 1.94 


13: 2.64x10-" 2.44 2.04 x 10- 1.98 
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Table 2: Numerical results for Nordsieck SGLMs of order p = q = 3. 


Type 2 method of order 3 [2] Nordsieck SGLM of order 3 


: en(Z) Pp en(Z) p 
4 4.07 x 10-6 6.58 x 1078 

5 5.66 x 10-7 2.85 8.66 x 1079 2.93 
6 7.43 x 1078 2.92 Lit 16-2 2.96 
7 9.50 x 10-9 2.96 1.40 x 1071° 2.99 


Table 3: Numerical results for Nordsieck SGLMs of order p = q = 4. 


Type 2 method of order 4 [2] Nordsieck SGLM of order 4 


/ en(Z) Pp en(Z) p 
4 6.22 x 1078 5.81 x 1079 

5 4.21 ~x 107-9 3.88 3.63 x 10-19 4.00 
6 2.74 10710 3.94 2.27 x 1071! 4.00 
7 175% i10-™ 3.97 1.42 x 10-1? 4.00 


4.2 Variable stepsize experiments 


We first investigate the potential for efficient implementation in a variable 
stepsize environment by the reliability of the estimation error for Problem 
1. The obtained results for the methods of order 2, 3, and 4 have plotted 
in Figures B, A, and , respectively. These figures confirm efficiency of the 
used estimation for the local truncation error. ‘To compare, we also present 
the results of numerical experiments of the LZ-stable Nordsieck SGLM given 
in [5]. In the implementation of considered SGLMs, we apply the same in- 
troduced implementation strategies, including the starting procedures, stage 
predictors, local error estimation, and the changing stepsize. In our numerical 
results, we use the following abbreviations: 


ns: the number of steps 

nrs: the number of rejected steps 

nfe: the number of function evaluations 
nJe: the number of Jacobian evaluations 
ge: the global error 

tol: given tolerance 


NSGLM2p: type 2 Nordsieck SGLM of order p 
NSGLM4p: type 4 Nordsieck SGLM of order p 
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Figure 3: Local errors and local error estimates versus x of the method of order 2 for 


problem 1 with ho = 10~5 and tol = 1078. 
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Figure 4: Local errors and local error estimates versus x of the method of order 3 for 


problem 1 with ho = 10~5 and tol = 1078. 
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Figure 5: Local errors and local error estimates versus x of the method of order 4 for 


problem 1 with ho = 10~5 and tol = 1078. 
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Table 4: Numerical results for Problem 2 solved by NSGLM2p and NSGLM4p of orders 
2,3, and 4 with ho = 107°. 


tol Method ge ns nrs nfe nJe 
10-° NSGLM22 3.61 x 10-° 169 26 1644 1256 

NSGLM42 2.02 x 10-4 1313 11 6693 4047 
10-? + NSGLM22 Lop i0-* 690 12 3718 2316 


NSGLM42 8.84 x 10~© 5997 14 24187 12167 


10-19 =NSGLM22 4.58 x 1078 3154 15 12847 6511 
NSGLM42 3.99 x 1077 27572 14 110499 55323 


10-° NSGLM23 2.95 x 107° 178 36 2339 1700 
NSGLM43 5.80 x 107° 426 17 3248 1922 
10-8 ~=NSGLM23 6.57 x 10-7 ATA 14 4047 2586 
NSGLM43 1.38 10-* 1341 23 8308 4219 
10-'° =NSGLM23 1.47 x 10-8 1450 18 8998 4597 
NSGLM43 5.59 x 10-9 26879 646 165153 82581 
10-° NSGLM24 9.16 x 10-° 337 25 3308 1864 
NSGLM44 9.66 x 10~° 313 14 3420 2116 
10-8 + NSGLM24 2.14 x 10~® 2040 50 16864 8508 
NSGLM44 2.19 x 10-® 787 14 6599 2399 


10-1 =NSGLM24 3.77 x 10-* 2696 141 22752 11408 
NSGLM44 4.59 x 10-8 1700 18 13934 7066 


Some numerical results for Problem 2 and 3 demonstrating the computational 
cost of NSGLM2p are given in Tables M and H and compared with those in 
NSGLM4p for p = 2,3,4. Also, in Figures § and [, we compare the accepted 
stepsizes of NSGLM23 with those in NSGLM24 and the accepted stepsizes of 
NSGLM43 with those in NSGLM44 through integration for Problem 2 and 
3, recpectively. The numerical results show that the proposed methods are 
capable in solving stiff problems and competitive with the existing methods. 


5 Conclusion 


We constructed type 2 Nordsieck SGLMs of orders 2, 3, and 4 with RKS 
property. Order 2 and 3 methods are L-stable and order 4 method is A- 
stable. These methods have been equipped to the variable stepsize using 
Nordsieck technique. The capability of the proposed methods in solving stiff 
problems with long interval of integration and badly scaled solution have 
been validated by some numerical experiments and comparisons. 
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Table 5: Numerical results for Problem 3 solved by NSGLM2p and NSGLM4p of orders 


2,3, and 4 with ho = 10~°®. 


tol Method ge ns nrs nfe nJe 
10-° NSGLM22 1.18 x 10~° 3336 492 24233 16579 
NSGLM42 8.95 x 107° 27540 A75 139942 83914 
10-8 NSGLM22 5.94 x 107-7 14721 486 77469 47057 
NSGLM42 3.94 x 10-® 126749 A427 508703 254353 
10-1 =NSGLM22 2.05 x 1078 67411 469 271519 135761 
NSGLM42 8.61 x 10-8 272719 396 1092456 546228 
10-° NSGLM23 itis 10-8 2790 A476 30994 21199 
NSGLM43 1.68 x 10~° 8243 499 61705 35482 
10-? + NSGLM23 1.46 x 107-7 8729 705 76184 74885 
NSGLM43 bol Ole 14444 450 155761 77884 
10-1 =NSGLM23 eae al teed 26879 646 165153 82581 
NSGLM43 Let <.10<* 80230 414 483861 241932 
10-° NSGLM24 Lex 105? 4070 642 42363 23519 
NSGLM44 2:11.x 10-3 4563 591 51425 30813 
10-8 =NSGLM24 4.86 x 107° 13440 1198 117655 59107 
NSGLM44 5.47 x 107° 10954 612 92614 46356 
10-1 =NSGLM24 8.86 x 10-7 24087 5512 236786 118394 
NSGLM44 £353¢10-% 26896 602 219980 109992 
10° 10” 


Figure 6: Accepted stepsizes versus x for Problem 2 with ho = 107% and tol = 1078: 
(a) NSGLM23 and NSGLM43, (b) NSGLM24 and NSGLM44. 
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—NSGLM23 —NSGLM24 
- - -NSGLM43 in [3] ~~ -NSGLM44 in [3] 


(a) (b) 


Figure 7: Accepted stepsizes versus « for Problem 3 with ho = 10~® and tol = 1078: 
(a) NSGLM23 and NSGLM43, (b) NSGLM24 and NSGLM44. 
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